library(tidyverse)
library(broom)
library(reshape2)
library(pbapply)
library(foreign)
library(stargazer)
library(lfe)
library(readstata13)
library(panelView)
library(data.table)
library(pbapply)
library(readxl)
library(haven)
library(labelled)
library(manifestoR)
library(lubridate)
library(scales)
library(ggplot2)
library(devtools)
library(psych)
library(xtable)
library(gmodels)
library(jackknifeR)
library(fixest)
library(scales)
library(gridExtra)

options(scipen = 999)
rm(list = ls())

##

############################################################
## Cross Country Analysis:
dat_comparative <- read_csv("dat_comparative_full.csv")

dat_comparative <- dat_comparative %>% 
  select(
    country_iso3c, election_date, election_id, election_monthyear, year, 
    turnout, polar_w_lag, disprop_lag, enp_lag, 
    east_dummy
  )

### ### ###
### Figure 2
dat_fig2 <- dat_comparative %>% select(election_date, year, country_iso3c, turnout, polar_w_lag ) %>% 
  group_by(country_iso3c) %>% 
  arrange(year, .by_group=T) %>%  ungroup() %>% 
  group_by(year) %>% 
  mutate(turnout_mean = mean(turnout),
         polar_mean = mean(polar_w_lag))  %>%  ungroup()
## Plot turnout over time in sample
ggplot(data=dat_fig2) +
  geom_smooth(aes(year, turnout), color = "grey35", linetype = "dashed") +
  geom_point(aes(x=year, y=turnout_mean)) +
  theme_bw() + ggplot2::theme(panel.grid.major = element_blank(), 
  panel.grid.minor = element_blank(), text = element_text(size = 15)) +
  ylab("Turnout") + xlab("Year")
ggsave("manuscript_figure_2.pdf", width = 16, height = 8, units = "cm")

### ### ###
### Table 1
m1_y <- felm(turnout ~ polar_w_lag + disprop_lag + enp_lag |  year | 0 | country_iso3c,
             data = dat_comparative)
m1_cy <- felm(turnout ~ polar_w_lag  + disprop_lag + enp_lag |  country_iso3c + year | 0 | country_iso3c,
              data = dat_comparative)

stargazer::stargazer( m1_y, m1_cy, 
                      type = "text",  omit = "factor",  omit.stat = c("rsq", "f", "ser"),
                      covariate.labels = c("Polarizationt t-1",
                                           "Disproportionality t-1",
                                           "ENP t-1"),
                      add.lines = list(c("Fixed Effects","Year", "Country+Year")),
                      out = c("manuscript_table_1.tex", "manuscript_table_1.html") )

### ### ###
### Table A7
m2_y <- felm(turnout ~ polar_w_lag  |  year | 0 | country_iso3c,
             data = dat_comparative)
m2_cy <- felm(turnout ~ polar_w_lag |  country_iso3c + year | 0 | country_iso3c,
              data = dat_comparative)

stargazer::stargazer( m2_y, m2_cy, 
                      type = "text",  omit = "factor",  omit.stat = c("rsq", "f", "ser"),
                      covariate.labels = c("Polarizationt t-1"),
                                           # ,
                                           # "Disproportionality t-1",
                                           # "ENP t-1"),
                      add.lines = list(c("Fixed Effects","Year", "Country+Year")),
                      out = c("manuscript_table_a7.tex", "manuscript_table_a7.html") )





### ### ### 
### Appendix Table A1
dat_descr <- dat_comparative %>% select(turnout, polar_w_lag, disprop_lag, enp_lag, year ) %>%  drop_na()
dat_descr_tab <- psych::describe(dat_descr)
dat_descr_tab$Variable <- row.names(dat_descr_tab)
dat_descr_tab <- as_tibble(dat_descr_tab) %>% relocate(Variable) %>% 
   dplyr::select(-vars, -trimmed, -skew, -kurtosis, -se, - mad)

print(xtable(dat_descr_tab), type="latex", include.rownames=FALSE, 
      file = "appendix_table_a1.tex" ) 


### ### ### 
### Appendix Table A2:
sample <- dat_fig2 %>% select(year, country_iso3c) %>% unique() %>%
  group_by(country_iso3c) %>% 
  summarise(Country = first(country_iso3c),
            min = min(year), 
            max = max(year)) %>%  select(-country_iso3c)
print(xtable(sample, digits = 0), type="latex", include.rownames=FALSE, file = "appendix_table_a2.tex" ) 


### ### ### 
### Table A8
jack_results_country1 <- tibble()
for(i in 1:length(dat_comparative$election_id) ){
  print(dat_comparative$election_id[i])
  country_id <- dat_comparative$country_iso3c[i]
  dat_jack <- dat_comparative %>% 
    filter( country_iso3c != country_id )
  
  m2_y_feols <- feols(turnout ~ polar_w_lag + disprop_lag + enp_lag | year,
                      cluster = ~country_iso3c, data = dat_jack)
  
  res_jack <- broom::tidy(m2_y_feols) %>% 
    mutate(excludes_country = country_id,
           n_obs_drop = 192 - m2_y_feols$nobs)
  
  jack_results_country1 <- bind_rows(jack_results_country1, res_jack) %>% distinct()
  }
jack_results_country1 <- jack_results_country1  %>% 
  filter(
    term == "polar_w_lag" ) %>% 
  select(-statistic) %>%   arrange(desc(p.value))

print(xtable(jack_results_country1, digits = 2), type="latex", include.rownames=FALSE, 
      file = "appendix_table_a8.tex" ) 

### ### ### 
## Table A9
jack_results_country2 <- tibble()
for(i in 1:length(dat_comparative$election_id)){
  print(dat_comparative$election_id[i])
  country_id <- dat_comparative$country_iso3c[i]
  dat_jack <- dat_comparative %>% 
    filter( country_iso3c != country_id  )
  
  m2_y_feols <- feols(turnout ~ polar_w_lag + disprop_lag + enp_lag | country_iso3c + year,
                      cluster = ~country_iso3c, data = dat_jack)
  
  res_jack <- broom::tidy(m2_y_feols) %>% 
    mutate(excludes_country = country_id,
           n_obs_drop = 192 - m2_y_feols$nobs)
  
  jack_results_country2 <- bind_rows(jack_results_country2, res_jack) %>% distinct()
}
jack_results_country2 <- jack_results_country2  %>% 
  filter(  term == "polar_w_lag"  ) %>% 
  select(-statistic) %>%  arrange(desc(p.value))

print(xtable(jack_results_country2, digits = 2), type="latex", include.rownames=FALSE, 
      file = "appendix_table_a9.tex" ) 


### ### ###
## Table a6
dat_cses <- read_xlsx("dat_CSES_PolarizationDatabase.xlsx") %>% 
  filter(Module == 4) 

dat_cses <- dat_cses %>% 
  group_by(ElectoralSystem) %>% 
  summarise( Polarization = mean(Polarization), 
             Cases = n(),
             # Country = any(Country)  
             Country = paste(unique(Country), collapse = ", ")
             )

print(xtable(dat_cses), type="latex", include.rownames=FALSE, 
      file = "appendix_table_a6.tex" ) 



##############################
## Diff-in-Diff Analyses
rm(list = ls())
dat_DE_county <- read_csv("dat_county_DEU_full_covs.csv")


### ### ### 
## Table 2:
m1 <- felm( turnout ~  afd_formed * afd_local | 0 | 0 | county_id_2018, data = dat_DE_county)
m2 <- felm( turnout ~  afd_formed * afd_local | state | 0 | county_id_2018, data = dat_DE_county)
m3 <- felm( turnout ~  afd_formed * afd_local | east | 0 | county_id_2018, data = dat_DE_county)
stargazer::stargazer( m1, m2, m3,
                      covariate.labels = c("AfD founded",
                                           "AfD candidate",
                                           "AfD founded*AfD candidate"),
                      type = "text", omit = "factor",  omit.stat = c("rsq", "f", "ser"), ci.level=0.95,
                      add.lines = list(c("FE?", "--", "State", "East")),
                      out = c("manuscript_table2.tex", "manuscript_table2.html"))



### ### ### 
## Appendix Table A12
dat_DE_county_east <- dat_DE_county %>%   filter(     east == 1   )
dat_DE_county_west <- dat_DE_county %>%   filter(    east == 0   )

m4 <- felm( turnout ~  afd_formed * afd_local | 0 | 0 | county_id_2018, data = dat_DE_county_west)
m5 <- felm( turnout ~  afd_formed * afd_local | 0 | 0 | county_id_2018, data = dat_DE_county_east)

stargazer::stargazer( m4, m5, 
                      covariate.labels = c("AfD founded",
                                           "AfD candidate",
                                           "AfD founded*AfD candidate"),
                       type = "text", omit = "factor",  omit.stat = c("rsq", "f", "ser"),
                       add.lines = list(c("Sample", "West", "East")),
                      out = c("appendix_table_a12.tex",  "appendix_table_a12.html"))


### ### ### 
## Appendix Table A3:
dat_descr <- dat_DE_county %>% select(turnout, afd_formed, afd_local, year )
dat_descr_tab <- psych::describe(dat_descr)   %>%   select(-vars, -trimmed, -skew, -kurtosis, -se, - mad)
print(xtable(dat_descr_tab), type="latex",  file = "appendix_table_a3.tex" ) 

dat_descr2 <- dat_DE_county %>% select(population_delta, gdppc_delta, emp_delta )
dat_descr_tab2 <- psych::describe(dat_descr2)   %>%   select(-vars, -trimmed, -skew, -kurtosis, -se, - mad)
print(xtable(dat_descr_tab2), type="latex",  file = "appendix_table_a3_covs.tex" ) 


### ### ### 
## Figure 3 (Parallel Trends):
dat_panel_full <- read_csv('dat_county_DEU_panel.csv')      
            
trend_1 <- ggplot(dat_panel_full, aes(x=year, y=turnout, color=as.factor(afd_local))) +
  scale_colour_hue(l=50) +  geom_smooth(method=lm, se=T ) +
  geom_vline(xintercept = 2013) +  labs(color = "AfD fielded\ncandidate") +
  theme_bw() + ggplot2::theme(panel.grid.major = element_blank(), 
                              panel.grid.minor = element_blank(), text = element_text(size = 15)) +
  theme(legend.position="bottom") + ylab("Turnout") + xlab("Year")
##
trend_2 <- ggplot(dat_panel_full, aes(x=year, y=turnout, color=as.factor(afd_local))) +
  scale_colour_hue(l=50) +  geom_smooth(method=loess) +
  geom_vline(xintercept = 2013) +   labs(color = "AfD fielded\ncandidate") +
  theme_bw() + ggplot2::theme(panel.grid.major = element_blank(), 
                              panel.grid.minor = element_blank(), text = element_text(size = 15)) +
  theme(legend.position="bottom") + ylab("Turnout") + xlab("Year")

trend_combined <- gridExtra::arrangeGrob(trend_1, trend_2, nrow = 1) 
ggsave("manuscript_figure_3.pdf", trend_combined,  width = 40, height = 20, units = "cm")


### ### ### 
## Table A13
dat_panel_full <- dat_panel_full %>%  mutate(afd_local = ifelse(year < 2013, 0, afd_local)) %>% 
  group_by(county_id_2018) %>%  arrange(year) %>%   mutate( turnout_delta = turnout - lag(turnout)) %>% ungroup()

m1 <- felm(afd_local  ~  turnout_delta  + year | 0 | 0 | county_id_2018, data = dat_panel_full)
m2 <- felm(afd_local  ~  turnout_delta  + year | state | 0 | county_id_2018,  data = dat_panel_full)

# Subset to main sample (pre-AfD)
dat_panel_subset <- dat_panel_full %>%  select(county_id_2018, year, turnout_delta)
dat_panel_subset <- left_join(dat_DE_county, dat_panel_subset) %>%  filter(year %in% c(2004:2013) )
m3 <- felm(afd_local  ~  turnout_delta  + year | 0 | 0 | county_id_2018, data = dat_panel_subset)

## Save 
stargazer::stargazer( m1, m2, m3,
                      covariate.labels = c("Delta Turnout", "Year"),
          type = "text", omit = "factor",  omit.stat = c("rsq", "f", "ser"),
                       add.lines = list(c("FE?", "--", "State","--" ), c("Sample", "Full", "Full" ,"2004-13" )),
                      out = c("appendix_table_a13.tex", "appendix_table_a13.html"))



### ### ###
## Table A14

# dat_DE_county<- read_csv("dat_county_DEU_full_covs.csv")

# Trends in pre-afd period
dat_DE_county_pre <- dat_DE_county %>% filter(afd_formed == 0) %>% ungroup() %>%  mutate(afd_local = as.numeric(afd_local))

t1 <- t.test( population_delta ~ afd_local , data=dat_DE_county_pre)
t1_df <- bind_cols("population_delta", as.numeric(t1$estimate[1]),as.numeric(t1$estimate[2]), as.numeric(t1$p.value) )

t2 <- t.test( gdppc_delta ~ afd_local , data=dat_DE_county_pre)
t2_df <- bind_cols("gdppc_delta", as.numeric(t2$estimate[1]), as.numeric(t2$estimate[2]), as.numeric(t2$p.value) )

t3 <- t.test( emp_delta ~ afd_local , data=dat_DE_county_pre)
t3_df <- bind_cols("emp_delta", as.numeric(t3$estimate[1]), as.numeric(t3$estimate[2]), as.numeric(t3$p.value) )

t_tests <- bind_rows( t1_df, t2_df, t3_df)
print(xtable(t_tests), type="latex", file = "appendix_table_a14.tex" ) 


### ### ###
## Table A3
dat_descr <- dat_DE_county %>% select(turnout, afd_formed, afd_local, year,
                                  population_delta, gdp_delta, emp_delta) %>%   mutate(year = as.numeric(year))

dat_descr_tab <- psych::describe(dat_descr) %>%   select(-vars, -trimmed, -skew, -kurtosis, -se, - mad)

print(xtable(dat_descr_tab), type="latex", file = "appendix_table_a3.tex" ) 


### ### ###
## Figure A1
m1 <- felm(afd_local  ~  population_delta | state | 0 | county_id_2018,
           data = dat_DE_county_pre) %>%   broom::tidy(conf.int = T)
m2 <- felm(afd_local  ~  gdppc_delta | state | 0 | county_id_2018,
           data = dat_DE_county_pre) %>%  broom::tidy(conf.int = T)
m3 <- felm(afd_local  ~  emp_delta | state | 0 | county_id_2018,
           data = dat_DE_county_pre) %>%  broom::tidy(conf.int = T)
results <- m1 %>% rbind(m2, m3) %>% 
    select(term,estimate,std.error,p.value,everything()) %>%   arrange( p.value)

dat_plot <- results %>% 
          mutate(  conf.low90 = estimate - qnorm(0.95)*std.error,
           conf.high90 = estimate + qnorm(0.95)*std.error,
           conf.low95 = estimate - qnorm(0.975) * std.error,
           conf.high95 = estimate + qnorm(0.975) * std.error)  %>% 
  mutate(term =  case_when(
    term == "population_delta" ~ "Change in population",
    term == "gdppc_delta" ~ "Change in GDP pc",
    term == "emp_delta" ~ "Change in employed constituents",
  ))

## PD
pd <- position_dodge(1)
## Plot
p1 <- ggplot(dat_plot[1,], aes(term, estimate)) +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_point(
    position = pd, shape = 21, fill  ='white') +
  geom_errorbar(aes(ymin = conf.low90, ymax = conf.high90),
                position = pd, width = 0, size = 1) +
  geom_errorbar(aes(ymin = conf.low95, ymax = conf.high95),
                                width = 0, size = 0.5) +
  geom_point( position = pd, shape = 21, fill  ='white') +
  theme(legend.position = 'bottom') +
  theme_bw() + ggplot2::theme(panel.grid.major = element_blank(), 
                              panel.grid.minor = element_blank(), text = element_text(size = 15)) +
  xlab('') +ylab('') +  theme(legend.position = 'bottom') 

p2 <- ggplot(dat_plot[2,], aes(term, estimate)) +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_point( position = pd,  shape = 21, fill  ='white') +
  geom_errorbar(aes(ymin = conf.low90, ymax = conf.high90),
                position = pd, width = 0, size = 1) +
  geom_errorbar(aes(ymin = conf.low95, ymax = conf.high95),
                width = 0, size = 0.5) +
  geom_point( position = pd,
    shape = 21,  fill  ='white') +
  theme(legend.position = 'bottom') +
  theme_bw() + ggplot2::theme(panel.grid.major = element_blank(), 
                              panel.grid.minor = element_blank(), text = element_text(size = 15)) +
  xlab('') +ylab('') +  theme(legend.position = 'bottom') 

p3 <- ggplot(dat_plot[3,], aes(term, estimate)) +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_point( position = pd, shape = 21, fill  ='white') +
  geom_errorbar(aes(ymin = conf.low90, ymax = conf.high90),
                position = pd, width = 0, size = 1) +
  geom_errorbar(aes(ymin = conf.low95, ymax = conf.high95),
                width = 0, size = 0.5) +
  geom_point( position = pd, shape = 21, fill  ='white') +
  theme(legend.position = 'bottom') +
  theme_bw() + ggplot2::theme(panel.grid.major = element_blank(), 
                              panel.grid.minor = element_blank(), text = element_text(size = 15)) +
  xlab('') +ylab('') + theme(legend.position = 'bottom') 
## Save
covar_plots <- arrangeGrob(p1, p2, p3, nrow = 1) 
ggsave("appendix_figure_a1.png", covar_plots,  width = 30, height = 15, units = "cm")


### ### ###
## Table a10
dat_county_polarization <- read_csv('dat_county_DEU_full_polarization.csv')

afd_local_prepost <- dat_county_polarization %>%
  filter(afd_local == 1) 
t1 <- t.test( polar_unit_alt ~ afd_formed , data=afd_local_prepost)
t1_df <- bind_cols("Polarization Index", as.numeric(t1$estimate[1]),as.numeric(t1$estimate[2]), as.numeric(t1$p.value) )

afd_local_post <- dat_county_polarization %>%
  filter(afd_formed == 1)

t2 <- t.test( polar_unit_alt ~ afd_local , data=afd_local_post)
t2
t2_df <- bind_cols("Polarization Index", as.numeric(t2$estimate[1]),as.numeric(t2$estimate[2]), as.numeric(t2$p.value) )

tab_a10 <- bind_rows(t1_df,t2_df)

print(xtable(tab_a10), type="latex", include.rownames=FALSE, 
      file = "appendix_table_a10.tex" ) 


### ### ###
## Table a11

dat_voting <- read_csv("dat_county_DEU_LocalNationalVoting.csv")


m1 <- felm( national_vs_extreme ~  county_vs_extreme | 0 | 0 | county_id_2018, data = dat_voting)
m2 <- felm( national_turnout ~  county_turnout | 0 | 0 | county_id_2018, data = dat_voting)

stargazer::stargazer( m1, m2,  
                      type = "text", omit = "factor",  omit.stat = c("rsq", "f", "ser"),
                      out = c("appendix_table_a11.tex", "appendix_table_a11.html"))



##################################################################
# Survey Analysis
rm(list = ls())

dat_survey <- read_csv("dat_CSES_survey.csv") 


### ### ###
# Table A15

m1 <- felm(turnout_pers ~ pos_pol_dalton + rile_congruence + 
             male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
           data = dat_survey)
m1_est <- broom::tidy(m1, conf.int=T)  %>%  mutate(model = "single", outcome = "Personal Vote", n = length(m1$residuals))   

m2 <- felm(turnout_pers ~ pos_sd + rile_congruence + 
             male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
           data = dat_survey)
m2_est <- broom::tidy(m2, conf.int=T)  %>%  mutate(model = "single",
                                                   outcome = "Personal Vote",
                                                   n = length(m2$residuals))   

m3 <- felm(turnout_pers ~ affect_high + rile_congruence + 
             male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
           data = dat_survey)
m3_est  <- broom::tidy(m3, conf.int=T)  %>% mutate(model = "single",
                                                   outcome = "Personal Vote",
                                                   n = length(m3$residuals))  

m4 <- felm(turnout_pers ~ affect_low + rile_congruence + 
             male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
           data = dat_survey)
m4_est <- broom::tidy(m4, conf.int=T)  %>% mutate(model = "single",
                                                  outcome = "Personal Vote",
                                                  n = length(m4$residuals))  


m1_nat <- felm(turnout_nat1 ~ pos_pol_dalton + rile_congruence + 
                 # party_id_str +
                 male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
               data = dat_survey)
m1nat_est <- broom::tidy(m1_nat, conf.int=T)  %>% mutate(model = "single",
                                                         outcome = "National Vote",
                                                         n = length(m1_nat$residuals)) 

m2_nat <- felm(turnout_nat1 ~ pos_sd + rile_congruence + 
                 # party_id_str +
                 male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
               data = dat_survey)
m2nat_est <- broom::tidy(m2_nat, conf.int=T)  %>% mutate(model = "single",
                                                         outcome = "National Vote",
                                                         n = length(m2_nat$residuals)) 

m3_nat <- felm(turnout_nat1 ~ affect_high + rile_congruence + 
                 # party_id_str +
                 male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
               data = dat_survey)

m3nat_est <- broom::tidy(m3_nat, conf.int=T)  %>% mutate(model = "single",
                                                         outcome = "National Vote",
                                                         n = length(m3_nat$residuals)) 

m4_nat <- felm(turnout_nat1 ~ affect_low + rile_congruence + 
                 # party_id_str +
                 male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
               data = dat_survey)
m4nat_est <- broom::tidy(m4_nat, conf.int=T)  %>% mutate(model = "single",
                                                         outcome = "National Vote",
                                                         n = length(m4_nat$residuals))  

stargazer::stargazer( m1, m2, m3, m4, m1_nat, m2_nat, m3_nat, m4_nat, 
                      type = "text", omit = "factor",  omit.stat = c("rsq", "f", "ser"),
                      ci.level=0.95, add.lines = list(c("Sociodem. Controls",  "Yes", "Yes", "Yes", "Yes", "Yes",  "Yes", "Yes", "Yes"),
                                                      c("C+Y FE?", "Yes", "Yes", "Yes","Yes", "Yes",  "Yes", "Yes", "Yes")),
                      out = c("appendix_table_a15.tex",
                              "appendix_table_a15.html"))

### ### ###
## Figure 4
ests <- bind_rows(m1_est, m2_est, m3_est, m4_est, 
                  m1nat_est, m2nat_est, m3nat_est, m4nat_est)
outcomes <- c("pos_pol_dalton" , "pos_sd" , "affect_high" , "affect_low")

ests <- ests %>% filter(term %in% outcomes) %>% 
  mutate(plot_lab = ifelse(model == "single", term, 'combined'))

## PD
pd <- position_dodge(1)
ests$term <- factor(ests$term, levels = c("pos_pol_dalton", "pos_sd", "affect_high", "affect_low"))
ests$term <- recode(ests$term, 
                    `pos_pol_dalton` = "Dalton Index",
                    `pos_sd` = "Spatial Deviation",
                    `affect_high` = "Positive Affect",
                    `affect_low` = "Negative Affect")
ests$plot_lab <- factor(ests$plot_lab, levels = c("pos_pol_dalton", "pos_sd", "affect_high", "affect_low", "combined"))
ests$plot_lab <- recode(ests$plot_lab, 
                        `pos_pol_dalton` = "Dalton Index",
                        `pos_sd` = "Spatial Deviation",
                        `affect_high` = "Positive Affect",
                        `affect_low` = "Negative Affect",
                        `combined` = "Combined Model")

ests$outcome <- factor(ests$outcome, levels = c("Personal Vote", "National Vote"))

## Plot
p1 <- ggplot(ests, aes(plot_lab, estimate, color = factor(term), shape = factor(term))) +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                position = pd,
                width = 0, size = 0.5) +
  geom_point(
    position = pd, fill  ='white', size = 3) +
  facet_wrap(~outcome, scales = 'free_y', ncol = 2) +
  theme_bw(base_size = 15) +
  xlab('Model Type (incl. Covariates and C+Y FEs)') +ylab('Coefficient (incl. 95% CIs)') +
  labs(color="Polarization Measure", shape = "Polarization Measure") +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90)) +
  theme_bw() + ggplot2::theme(panel.grid.major = element_blank(), 
                              panel.grid.minor = element_blank(), text = element_text(size = 15)) +
  theme(legend.position = 'bottom') 

p1
ggsave('manuscript_figure_4.pdf', width = 28, height = 15, units = "cm")




### ### ###
## Table A16 

m4_affect_extreme <- felm(affect_low ~   rile_extreme +
                            male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
                          data = dat_survey)

m4_affect_leftright <- felm(affect_low ~   rile_left +  rile_right +
                              male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
                            data = dat_survey)

m4_affect_knowledge <- felm(affect_low ~   know_high + know_low +
                              male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
                            data = dat_survey)


stargazer::stargazer( m4_affect_extreme, m4_affect_leftright, m4_affect_knowledge, 
                      keep = c("rile_extreme", "rile_left", "rile_right", "know_high", "know_low"),
                      type = "text", omit = "factor",  omit.stat = c("rsq", "f", "ser"),
                      add.lines = list(c("Sociodem. Controls",  "Yes", "Yes", "Yes", "Yes"),
                                       c("C+Y FE?", "Yes", "Yes", "Yes","Yes")),
                      out = c("appendix_table_a16.tex",  "appendix_table_a16.html"))

### ### ###
## Tables A17
m4_close_nat <- felm(turnout_nat1 ~ affect_low +  closeness_votes +
                       male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
                     data = dat_survey)

m4_close_pers <- felm(turnout_pers ~ affect_low +  closeness_votes +
                        male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
                      data = dat_survey)

m4_close_affect <- felm(affect_low ~   closeness_votes +
                          male + age + income + education  | ccode_un + eyear | 0 | ccode_un,
                        data = dat_survey)


stargazer::stargazer( m4_close_affect, m4_close_nat, m4_close_pers, 
                      keep = c("affect_low", "closeness_votes"),
                      type = "text", omit = "factor",  omit.stat = c("rsq", "f", "ser"),
                      add.lines = list(c("Individual Controls",  "Yes", "Yes", "Yes"),
                                       c("C+Y FE?", "Yes", "Yes", "Yes")),
                      out = c("appendix_table_a17.tex",  "appendix_table_a17.html"))


### ### ###
## Table a4
dat_descr <- dat_survey %>% select(turnout_nat1, turnout_pers,
                            pos_pol_dalton, pos_sd, affect_high, affect_low, 
                            closeness_votes,rile_congruence,
                            rile_extreme, rile_left, rile_right,
                            know_high, know_low, male, age, income, education, eyear )

dat_descr_tab <- psych::describe(dat_descr)   %>%   select(-vars, -trimmed, -skew, -kurtosis, -se, - mad)

print(xtable(dat_descr_tab), type="latex",
      file = "appendix_table_a4.tex" ) 


### ### ###
## Table a5

countries <- dat_survey %>% dplyr::select(cname, eyear) %>%   group_by(cname) %>% 
  summarise(Min = min(eyear), Max = max(eyear))
print(xtable(countries, digits = 0), type="latex", include.rownames=FALSE, 
      file = "appendix_table_a5.tex" ) 





